### Please also see README for required R packages
packages <- c(
  "tidyverse",
  "firatheme",
  "ggplot2",
  "Hmisc",
  "broom",
  "MASS",
  "effects",
  "cjoint",
  "scales",
  "stringr",
  "forcats",
  "patchwork",
  "gridExtra",
  "cowplot"
)

lapply(packages, library, character.only = TRUE)

### load data
load("data/fig2/fig2.RData")

load("data/fig3/embeddings_w1.RData")
load("data/fig3/embeddings_w2.RData")
load("data/fig3/embeddings_w3.RData")
load("data/fig3/embeddings_pooled.RData")

load("data/fig3/preprocessed_w1.RData")
load("data/fig3/preprocessed_w2.RData")
load("data/fig3/preprocessed_w3.RData")
load("data/fig3/preprocessed_pooled.RData")

load("data/fig3/rawtext_w1.RData")
load("data/fig3/rawtext_w2.RData")
load("data/fig3/rawtext_w3.RData")
load("data/fig3/rawtext_pooled.RData")

load("data/fig3/distance_w1.RData")
load("data/fig3/distance_w2.RData")
load("data/fig3/distance_w3.RData")
load("data/fig3/distance_pooled.RData")

load("data/fig4/class_w1.RData")
load("data/fig4/class_w2.RData")
load("data/fig4/class_w3.RData")
load("data/fig4/class_pooled.RData")

load("data/fig4/education_w1.RData")
load("data/fig4/education_w2.RData")
load("data/fig4/education_w3.RData")
load("data/fig4/education_pooled.RData")

load("data/fig5/correct_answer.RData")
load("data/fig5/gender.RData")
load("data/fig5/party.RData")

load("data/fig6/conjoint_low_information.RData")
load("data/fig7/conjoint_high_information.RData")

load("data/fig8/conjoint_low_information.RData")
load("data/fig8/conjoint_high_information.RData")


### Figure 2 - upper left panel
boxplot(W~condition, data=fig2, notch=F, ylab = "",
        col=(c("gray70","gray90")),
        names = c("high", "low"), main="Number of words", xlab="Sophistication")

### Figure 2 - upper center panel 
boxplot(C~condition, data=fig2, notch=F, ylab = "",
        col=(c("gray70","gray90")),
        main="Number of characters", xlab="Sophistication",
        names = c("high", "low"), ylim=c(200,460))

### Figure 2 - upper right panel 
boxplot(meanWordChars~condition, data=fig2, notch=F, ylab = "",
        col=(c("gray70","gray90")),
        main="Avergage characters \nper word", xlab="Sophistication",
        names = c("high", "low"), ylim=c(5,7.5))

### Figure 2 - lower left panel 
boxplot(meanSentenceChars~condition, data=fig2, ylab = "",
        col=(c("gray70","gray90")),
        main="Avergage characters \nper sentence", xlab="Sophistication",
        names = c("high", "low"))

### Figure 2 - lower center panel 
boxplot(meanSentenceLength~condition, data=fig2,ylab = "",
        col=(c("gray70","gray90")),
        main="Average number of \nwords per sentence", xlab="Sophistication",
        names = c("high", "low"))

### Figure 2 - lower right panel 
boxplot(per_leipzig~condition, data=fig2, ylab = "",
        col=(c("gray70","gray90")),
        main="Percentage Leipzig \nCorpus Top 1000", xlab="Sophistication",
        names = c("high", "low"), ylim=c(55,70))

### Figure 3 - left panel - embeddings
embeddings_model_w1 <- lm(cosine_scores ~ treatment + issue + position, data = embeddings_w1)
embeddings_model_w2 <- lm(cosine_scores ~ treatment + issue + position, data = embeddings_w2)
embeddings_model_w3 <- lm(cosine_scores ~ treatment + issue + position, data = embeddings_w3)
embeddings_model_pooled <- lm(cosine_scores ~ treatment + issue + position + fe, data = embeddings_pooled)

### Figure 3 - left panel - pre-processed
preprocessed_model_w1 <- lm(cosine ~ treatment + issue + position, data = preprocessed_w1)
preprocessed_model_w2 <- lm(cosine ~ treatment + issue + position, data = preprocessed_w2)
preprocessed_model_w3 <- lm(cosine ~ treatment + issue + position, data = preprocessed_w3)
preprocessed_model_pooled <- lm(cosine ~ treatment + issue + position + fe, data = preprocessed_pooled)

### Figure 3 - left panel - raw text
rawtext_model_w1 <- lm(cosine ~ treatment + issue + position, data = rawtext_w1)
rawtext_model_w2 <- lm(cosine ~ treatment + issue + position, data = rawtext_w2)
rawtext_model_w3 <- lm(cosine ~ treatment + issue + position, data = rawtext_w3)
rawtext_model_pooled <- lm(cosine ~ treatment + issue + position + fe, data = rawtext_pooled)

# Function to extract coefficient and confidence intervals
extract_ci <- function(model, variable) {
  ci <- confint(model, parm = variable)  # Get confidence intervals
  coef <- coef(summary(model))[variable, "Estimate"]  # Get coefficient
  data.frame(
    ate = coef,
    lower = ci[1],
    upper = ci[2]
  )
}

# Variable of interest (treatment effect)
treatment_var <- "treatment"

# Extract estimates from each model
df_list <- list(
  extract_ci(embeddings_model_w1, treatment_var) %>% mutate(wave = "wave1", type2 = "3 embeddings", type = "wave estimates"),
  extract_ci(embeddings_model_w2, treatment_var) %>% mutate(wave = "wave2", type2 = "3 embeddings", type = "wave estimates"),
  extract_ci(embeddings_model_w3, treatment_var) %>% mutate(wave = "wave3", type2 = "3 embeddings", type = "wave estimates"),
  extract_ci(embeddings_model_pooled, treatment_var) %>% mutate(wave = "pooled", type2 = "3 embeddings", type = "pooled estimate"),
  
  extract_ci(preprocessed_model_w1, treatment_var) %>% mutate(wave = "wave1", type2 = "2 pre-processed", type = "wave estimates"),
  extract_ci(preprocessed_model_w2, treatment_var) %>% mutate(wave = "wave2", type2 = "2 pre-processed", type = "wave estimates"),
  extract_ci(preprocessed_model_w3, treatment_var) %>% mutate(wave = "wave3", type2 = "2 pre-processed", type = "wave estimates"),
  extract_ci(preprocessed_model_pooled, treatment_var) %>% mutate(wave = "pooled", type2 = "2 pre-processed", type = "pooled estimate"),
  
  extract_ci(rawtext_model_w1, treatment_var) %>% mutate(wave = "wave1", type2 = "1 raw text", type = "wave estimates"),
  extract_ci(rawtext_model_w2, treatment_var) %>% mutate(wave = "wave2", type2 = "1 raw text", type = "wave estimates"),
  extract_ci(rawtext_model_w3, treatment_var) %>% mutate(wave = "wave3", type2 = "1 raw text", type = "wave estimates"),
  extract_ci(rawtext_model_pooled, treatment_var) %>% mutate(wave = "pooled", type2 = "1 raw text", type = "pooled estimate")
)

# Combine into a single data frame
fig3_data <- bind_rows(df_list)

# Plot left panel of Figure 3
ggplot(fig3_data, aes(wave, ate, group = type2)) +
  geom_point(aes(color = type2), stroke = 0, position = position_dodge(width = 0.5)) +
  geom_linerange(aes(ymin = lower, ymax = upper, color = type2), position = position_dodge(width = 0.5)) +
  geom_hline(yintercept = 0, linetype = 1, color = "black", alpha = 0.4) +
  facet_wrap(~ type, scales = "free_x") +
  theme_fira() +
  ylab("Cosine similarity \n(scale 0-1)") +
  ylim(-0.1, 0.012) +
  scale_color_manual(values = c("black", "gray40", "gray80")) +  # Set black and gray colors
  theme(
    strip.background = element_blank(),
    axis.title.x = element_blank(),
    panel.grid.minor = element_blank(),
    text = element_text(size = 8),
    legend.position = "bottom",
    legend.box.background = element_rect(fill = "white"),
    legend.key.size = unit(1, "lines"),
    legend.direction = "horizontal",
    legend.title = element_blank(),
    strip.text.x = element_text(size = 9, colour = "black")
  )


### Figure 3 - right panel 
distance_model_w1 <- lm(outcome2_distance ~ treatment + issue + position, data = distance_w1)
distance_model_w2 <- lm(outcome2_distance ~ treatment + issue + position, data = distance_w2)
distance_model_w3 <- lm(outcome2_distance ~ treatment + issue + position, data = distance_w3)
distance_model_pooled <- lm(outcome2_distance ~ treatment + issue + position + fe, data = distance_pooled)

# Run Code: Function to extract coefficient and confidence intervals: Line 120-129 

# Extract estimates from each model
df_list <- list(
  extract_ci(distance_model_w1, treatment_var) %>% mutate(wave = "wave1", type = "wave estimates"),
  extract_ci(distance_model_w2, treatment_var) %>% mutate(wave = "wave2", type = "wave estimates"),
  extract_ci(distance_model_w3, treatment_var) %>% mutate(wave = "wave3", type = "wave estimates"),
  extract_ci(distance_model_pooled, treatment_var) %>% mutate(wave = "pooled", type = "pooled estimate")
)

# Combine into a single data frame
fig_distance <- bind_rows(df_list)

# Plot
ggplot(fig_distance, aes(wave, ate)) +
  geom_point(stroke = 0, size = 2) +
  geom_linerange(aes(ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0, linetype = 1, color = "black", alpha = 0.4) +
  facet_wrap(~ type, scales = "free_x") +
  theme_fira() +
  ylab("Distance to correct answer \n(scale 0-5)") +
  ylim(-0.15, 0.3) +
  theme(
    strip.background = element_blank(),
    axis.title.x = element_blank(),
    panel.grid.minor = element_blank(),
    text = element_text(size = 8),
    legend.position.inside = c(0.6, 0.06),
    legend.box.background = element_rect(fill = "gray"),
    legend.key.size = unit(1, "lines"),
    legend.direction = "horizontal",
    legend.title = element_blank(),
    strip.text.x = element_text(size = 9, colour = "black")
  )

### Figure 4 - left panel
class_model_w1 <- polr(class_order ~ treatment + issue + position, data = class_w1)
class_model_w2 <- polr(class_order ~ treatment + issue + position, data = class_w2)
class_model_w3 <- polr(class_order ~ treatment + issue + position, data = class_w3)
class_model_pooled <- polr(class_order ~ treatment + issue + position + fe, data = class_pooled)


# Function to extract predicted probabilities using Effect()
extract_class_effects <- function(model, wave_label) {
  effects <- Effect(focal.predictors = "treatment", model)
  
  # Convert effect data to a dataframe
  effect_data <- as.data.frame(effects)
  
  # Ensure treatment is correctly assigned as a factor
  effect_data$treatment <- as.factor(effect_data$treatment)  
  levels(effect_data$treatment) <- c("simple", "sophisticated")  # Rename levels
  
  # Extract relevant probabilities (dropping Mittelschicht)
  extracted_data <- data.frame(
    prob = c(effect_data[["prob.Unterschicht"]], effect_data[["prob.Oberschicht"]]),  # Lower & Upper class
    lower = c(effect_data[["L.prob.Unterschicht"]], effect_data[["L.prob.Oberschicht"]]), 
    upper = c(effect_data[["U.prob.Unterschicht"]], effect_data[["U.prob.Oberschicht"]]),
    group = rep(levels(effect_data$treatment), 2),  # "simple" and "sophisticated"
    category = rep(c("lower class", "upper class"), each = 2),  # Assign class categories
    wave = wave_label,
    type = ifelse(wave_label == "pooled", "pooled estimates", "wave estimates")
  )
  
  return(extracted_data)
}

# Extract effects for each model
df_list <- list(
  extract_class_effects(class_model_w1, "wave1"),
  extract_class_effects(class_model_w2, "wave2"),
  extract_class_effects(class_model_w3, "wave3"),
  extract_class_effects(class_model_pooled, "pooled")
)

# Combine all into a single data frame
fig3_data <- bind_rows(df_list)

# Ensure category ordering
fig3_data$category <- factor(fig3_data$category, levels = c("upper class", "lower class"))

# Assign different colors for simple and sophisticated groups
color_mapping <- c("simple" = "gray60", "sophisticated" = "black")

# Plot
ggplot(fig3_data, aes(wave, prob, color = group)) +
  geom_point(stroke = 0, position = position_dodge(width = 0.5)) +
  geom_linerange(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.5), alpha = 0.8) +
  scale_y_continuous(expand = c(0.05, 0.05)) +
  facet_grid(category ~ type, scales = "free") + 
  theme_fira() +
  ylab("Probability") +
  theme(
    strip.background = element_blank(),
    axis.title.x = element_blank(),
    panel.grid.minor = element_blank(),
    text = element_text(size = 8),
    legend.position = "bottom",
    legend.background = element_rect(linewidth = 0.5), 
    legend.key.size = unit(1, "lines"),
    legend.direction = "horizontal",
    legend.title = element_blank(),
    legend.text = element_text(size = 5),
    strip.text.x = element_text(size = 9, colour = "black")
  ) +
  scale_color_manual(values = color_mapping)  # Assign correct colors for each group

### Figure 4 - right panel
education_model_w1 <- polr(edu_order ~ treatment + issue + position, data = education_w1)
education_model_w2 <- polr(edu_order ~ treatment + issue + position, data = education_w2)
education_model_w3 <- polr(edu_order ~ treatment + issue + position, data = education_w3)
education_model_pooled <- polr(edu_order ~ treatment + issue + position + fe, data = education_pooled)

# Function to extract predicted probabilities using Effect()
extract_education_effects <- function(model, wave_label) {
  effects <- Effect(focal.predictors = "treatment", model)
  
  # Convert effect data to a dataframe
  effect_data <- as.data.frame(effects)
  
  # Ensure treatment is correctly assigned as a factor
  effect_data$treatment <- as.factor(effect_data$treatment)  
  levels(effect_data$treatment) <- c("simple", "sophisticated")  # Rename levels
  
  # Extract relevant probabilities (dropping "mittlere Bildung")
  extracted_data <- data.frame(
    prob = c(effect_data[["prob.niedrige.Bildung"]], effect_data[["prob.hohe.Bildung"]]),  # Low & High education
    lower = c(effect_data[["L.prob.niedrige.Bildung"]], effect_data[["L.prob.hohe.Bildung"]]), 
    upper = c(effect_data[["U.prob.niedrige.Bildung"]], effect_data[["U.prob.hohe.Bildung"]]),
    group = rep(levels(effect_data$treatment), 2),  # "simple" and "sophisticated"
    category = rep(c("low education", "high education"), each = 2),  # Assign education categories
    wave = wave_label,
    type = ifelse(wave_label == "pooled", "pooled estimates", "wave estimates")
  )
  
  return(extracted_data)
}

# Extract effects for each model
df_list <- list(
  extract_education_effects(education_model_w1, "wave1"),
  extract_education_effects(education_model_w2, "wave2"),
  extract_education_effects(education_model_w3, "wave3"),
  extract_education_effects(education_model_pooled, "pooled")
)

# Combine all into a single data frame
fig4_data <- bind_rows(df_list)

# Ensure category ordering
fig4_data$category <- factor(fig4_data$category, levels = c("high education", "low education"))

# Assign different colors for simple and sophisticated groups
color_mapping <- c("simple" = "gray60", "sophisticated" = "black")

# Plot
ggplot(fig4_data, aes(wave, prob, color = group)) +
  geom_point(stroke = 0, position = position_dodge(width = 0.5)) +
  geom_linerange(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.5), alpha = 0.8) +
  scale_y_continuous(expand = c(0.05, 0.05)) +
  facet_grid(category ~ type, scales = "free") + 
  theme_fira() +
  ylab("Probability") +
  theme(
    strip.background = element_blank(),
    axis.title.x = element_blank(),
    panel.grid.minor = element_blank(),
    text = element_text(size = 8),
    legend.position = "bottom",
    legend.background = element_rect(size = 0.5), 
    legend.key.size = unit(1, "lines"),
    legend.direction = "horizontal",
    legend.title = element_blank(),
    legend.text = element_text(size = 5),
    strip.text.x = element_text(size = 9, colour = "black")
  ) +
  scale_color_manual(values = color_mapping)  # Assign correct colors for each group

### Figure 5 
model_a <- glm(profile_gender_correct ~ treatment + issue + position + fe, data = data_correct_answer)
model_b <- glm(profile_party_correct ~ treatment + issue + position + fe, data = data_correct_answer)
model_c <- polr(profil_party_0_GROUP ~ treatment + issue + position + fe, data = data_party)
model_d <- polr(profil_gender_0_GROUP ~ treatment + issue + position + fe, data = data_gender)

extract_glm_effects <- function(model, category_label) {
  effects <- Effect(focal.predictors = "treatment", model)
  effect_data <- as.data.frame(effects)
  
  # Ensure treatment is a factor and rename levels
  effect_data$treatment <- factor(effect_data$treatment, levels = c("simple", "sophisticated"))
  
  # Extract relevant probabilities
  extracted_data <- data.frame(
    prob = effect_data$fit,
    lower = effect_data$lower,
    upper = effect_data$upper,
    group = effect_data$treatment,
    category = category_label
  )
  
  return(extracted_data)
}

# Extract GLM effects
glm_data <- bind_rows(
  extract_glm_effects(model_a, "correct gender"),
  extract_glm_effects(model_b, "correct party")
)

extract_polr_party_effects <- function(model) {
  effects <- Effect(focal.predictors = "treatment", model)
  effect_data <- as.data.frame(effects)
  
  # Ensure treatment is a factor and rename levels
  effect_data$treatment <- factor(effect_data$treatment, levels = c("simple", "sophisticated"))
  
  # Define correct column names **(matching your output)**
  categories <- c("No Party", "AfD", "Grüne", "CDU/CSU", "Linke", "FDP", "SPD")
  raw_labels <- c("prob.X.99",  
                  "prob.AfD..Die.Alternative.für.Deutschland.", 
                  "prob.Bündnis.90.Die.Grünen", 
                  "prob.CDU.CSU..Christlich.Demokratische.Union.Deutschlands.", 
                  "prob.Die.Linke", 
                  "prob.FDP..Freie.Demokratische.Partei.", 
                  "prob.SPD..Sozialdemokratische.Partei.Deutschlands.")
  
  lower_labels <- sub("prob", "L.prob", raw_labels)
  upper_labels <- sub("prob", "U.prob", raw_labels)
  
  # Extract relevant probabilities **(only the relevant values)**
  extracted_data <- data.frame(
    prob = unlist(lapply(raw_labels, function(col) effect_data[[col]])),
    lower = unlist(lapply(lower_labels, function(col) effect_data[[col]])),
    upper = unlist(lapply(upper_labels, function(col) effect_data[[col]])),
    group = rep(levels(effect_data$treatment), length(categories)),  
    category = rep(categories, each = 2)  # Two values per category (one per treatment)
  ) 
  
  return(extracted_data)
}


# Extract POLR (party) effects with correct labels
polr_party_data <- extract_polr_party_effects(model_c)
polr_party_data <- polr_party_data |>
  filter(category != "No Party")

extract_polr_gender_effects <- function(model) {
  effects <- Effect(focal.predictors = "treatment", model)
  effect_data <- as.data.frame(effects)
  
  # Ensure treatment is a factor and rename levels
  effect_data$treatment <- factor(effect_data$treatment, levels = c("simple", "sophisticated"))
  
  # Define gender categories **(fixing to "weiblich" & "männlich")**
  categories <- c("male", "female")
  
  # Extract relevant probabilities **(ensuring correct labels)**
  extracted_data <- data.frame(
    prob = c(effect_data[["prob.männlich"]], effect_data[["prob.weiblich"]]),
    lower = c(effect_data[["L.prob.männlich"]], effect_data[["L.prob.weiblich"]]),
    upper = c(effect_data[["U.prob.männlich"]], effect_data[["U.prob.weiblich"]]),
    group = rep(levels(effect_data$treatment), 2),
    category = rep(categories, each = 2)
  )
  
  return(extracted_data)
}

# Extract POLR (gender) effects
polr_gender_data <- extract_polr_gender_effects(model_d)

# Add id variable to each dataset
glm_data <- glm_data |> mutate(id = "correct answer")
polr_party_data <- polr_party_data |> mutate(id = "party")
polr_gender_data <- polr_gender_data |> mutate(id = "gender")

# Combine all datasets
combined_data <- bind_rows(glm_data, polr_party_data, polr_gender_data)

# Create the faceted plot
ggplot(combined_data, aes(category, prob, color = group)) +
  geom_point(stroke = 0, position = position_dodge(width = 0.5)) +
  geom_linerange(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.5), alpha = 0.8) +
  scale_y_continuous(expand = c(0.05, 0.05)) +
  facet_wrap(~id, scales = "free") +  # Facet by dataset type
  theme_fira() +
  ylab("Probability") +
  theme(
    strip.background = element_blank(),
    axis.title.x = element_blank(),
    panel.grid.minor = element_blank(),
    text = element_text(size = 8),
    legend.position = "bottom",
    legend.background = element_rect(size = 0.5), 
    legend.key.size = unit(1, "lines"),
    legend.direction = "horizontal",
    legend.title = element_blank(),
    legend.text = element_text(size = 8),
    strip.text.x = element_text(size = 9, colour = "black")
  ) +
  scale_color_manual(values = c("gray60", "black"))


### Figure 6 and 7 

# Function to extract AMCE results
extract_amce <- function(model, dataset_label) {
  summary_data <- summary(model)
  
  # Extract AMCE estimates and baselines
  amce_data <- summary_data$amce |> dplyr::select(Attribute, Level, Estimate, `Std. Err`, `z value`)
  baselines_data <- summary_data$baselines_amce |> 
    mutate(Estimate = NA, `Std. Err` = NA, `z value` = NA)
  
  # Combine and compute confidence intervals
  bind_rows(amce_data, baselines_data) |>
    mutate(
      upper = Estimate + 1.96 * `Std. Err`,
      lower = Estimate - 1.96 * `Std. Err`,
      dataset = dataset_label
    ) |>
    replace_na(list(Estimate = 0, upper = 0, lower = 0))
}

# Low Information AMCE
amce2 <- cjoint::amce(selected ~ mp_sex + mp_age + mp_party + mp_position + mp_sophistication, 
                      cjointdata2, design = MPdesign2, respondent.id = "respondent")
plotdata_amce2 <- extract_amce(amce2, "Low Information")

# High Information AMCE
amce1 <- cjoint::amce(selected ~ mp_sex + mp_age + mp_party + mp_motivation + mp_days +  
                        mp_occupation + mp_experience + mp_position + mp_sophistication, 
                      cjointdata1, design = MPdesign, respondent.id = "respondent")
plotdata_amce1 <- extract_amce(amce1, "High Information")

# Combine All Data
combined_amce <- bind_rows(plotdata_amce1, plotdata_amce2)

# Define Order and Baselines More Efficiently
combined_amce <- combined_amce |>
  mutate(
    plotorder = factor(Attribute, levels = c("mp_age", "mp_party", "mp_sex", "mp_position", "mp_sophistication",
                                             "mp_motivation", "mp_days", "mp_experience", "mp_occupation"),
                       labels = c("age", "party", "sex", "position", "sophistication", 
                                  "motivation", "days", "experience", "occupation")),
    
    # Define Order Efficiently
    order = case_when(
      Level %in% c("26", "AfD", "against", "male", "simple", "to serve the party", "1 day", "new", "farmer") ~ 1,
      Level %in% c("30", "Die Linke", "neutral", "female", "sophisticated", "to impact personally", "2 days", "incumbent", "nurse") ~ 2,
      Level %in% c("47", "B90/Grüne", "for", "to represent ordinary people", "3 days", "nurturer") ~ 3,
      Level %in% c("57", "FDP", "lawyer") ~ 4,
      Level %in% c("64", "CDU/CSU", "political scientist") ~ 5,
      Level %in% c("79", "SPD", "teacher") ~ 6,
      TRUE ~ NA_real_
    ),
    
    # Define Baselines More Efficiently
    baseline = if_else(Level %in% c("AfD", "against", "male", "simple", "26", "new", "1 day", "farmer", "to serve the party"), "1", "0"),
    
    # Shorten Labels for Better Readability
    Level = dplyr::recode(Level, 
                   "to serve the party" = "serve party", 
                   "to impact personally" = "impact personally", 
                   "to represent ordinary people" = "represent people")
  )

plot_amce <- function(data, title_text) {
  ggplot(data, aes(x = reorder(Level, -order), y = Estimate, color = baseline)) +
    geom_hline(yintercept = 0, linetype = 2, color = "gray70") +
    geom_pointrange(aes(ymin = lower, ymax = upper), size = 0.7, fatten = 0.5, position = position_dodge(width = 0.8)) +
    coord_flip() +
    scale_y_continuous(limits = c(-0.15, 0.3), breaks = c(-0.1, 0, 0.1, 0.2)) +
    facet_wrap(~plotorder, scales = "free", ncol = 3) +
    scale_colour_manual(values = c("gray30", "gray80")) +
    theme_fira() +
    labs(title = title_text, x = "", y = "") +
    theme(legend.position = "none",
          axis.text.y = element_text(size = 14, color = "black"),
          axis.text.x = element_text(size = 14, color = "black"),
          strip.text = element_text(size = 14, color = "black"))
}

# Generate Plots - Figure 6 and 7 
plot_amce(filter(combined_amce, dataset == "Low Information"), "Low-Information Environment")
plot_amce(filter(combined_amce, dataset == "High Information"), "High-Information Environment")


# Figure 8 
# Function to extract AMCE results and filter only `mp_position`
extract_mp_position <- function(model, dataset_label) {
  summary_data <- summary(model)
  
  # Extract AMCE estimates and baselines, keeping only `mp_position`
  amce_data <- summary_data$amce |> 
    dplyr::select(Attribute, Level, Estimate, `Std. Err`, `z value`) |> 
    filter(Attribute == "mp_position")
  
  baselines_data <- summary_data$baselines_amce |> 
    mutate(Estimate = NA, `Std. Err` = NA, `z value` = NA) |> 
    filter(Attribute == "mp_position")
  
  # Combine and compute confidence intervals
  bind_rows(amce_data, baselines_data) |>
    mutate(
      upper = Estimate + 1.96 * `Std. Err`,
      lower = Estimate - 1.96 * `Std. Err`,
      dataset = dataset_label
    ) |>
    replace_na(list(Estimate = 0, upper = 0, lower = 0))
}

# Low Information (Separate by Sophistication)
cjointdata2_sophisticated <-  subset(cjointdata2, mp_sophistication == "sophisticated")
cjointdata2_simple        <-  subset(cjointdata2, mp_sophistication == "simple")

amce2_sophisticated <- cjoint::amce(selected ~ mp_sex + mp_age + mp_party + mp_position, 
                                    cjointdata2_sophisticated, design = MPdesign2, respondent.id = "respondent")

amce2_simple <- cjoint::amce(selected ~ mp_sex + mp_age + mp_party + mp_position, 
                             cjointdata2_simple, design = MPdesign2, respondent.id = "respondent")

plotdata_amce2_soph <- extract_mp_position(amce2_sophisticated, "Low-Info Sophisticated")
plotdata_amce2_simple <- extract_mp_position(amce2_simple, "Low-Info Simple")

# High Information (Separate by Sophistication)
cjointdata1_sophisticated <-  subset(cjointdata1, mp_sophistication == "sophisticated")
cjointdata1_simple        <-  subset(cjointdata1, mp_sophistication == "simple")

amce1_sophisticated <- cjoint::amce(selected ~ mp_sex + mp_age + mp_party + mp_motivation + mp_days +  
                                      mp_occupation + mp_experience + mp_position, 
                                    cjointdata1_sophisticated, design = MPdesign, respondent.id = "respondent")

amce1_simple <- cjoint::amce(selected ~ mp_sex + mp_age + mp_party + mp_motivation + mp_days +  
                               mp_occupation + mp_experience + mp_position, 
                             cjointdata1_simple, design = MPdesign, respondent.id = "respondent")

plotdata_amce1_soph <- extract_mp_position(amce1_sophisticated, "High-Info Sophisticated")
plotdata_amce1_simple <- extract_mp_position(amce1_simple, "High-Info Simple")

# Combine Data
mp_position_data <- bind_rows(plotdata_amce1_soph, plotdata_amce1_simple, 
                              plotdata_amce2_soph, plotdata_amce2_simple)

mp_position_data <- mp_position_data |>
  mutate(
    plotorder = factor(Attribute, 
                       levels = c("mp_position"),  # Only keep mp_position
                       labels = c("position")),  # Label it as 'position'
    
    order = case_when(
      Level == "against" ~ 1,
      Level == "neutral" ~ 2,
      Level == "for" ~ 3,
      TRUE ~ NA_real_
    ),
    
    baseline = if_else(Level == "against", "1", "0")
  )

plot_mp_position <- function(data, title_text) {
  ggplot(data, aes(x = reorder(Level, -order), y = Estimate, color = dataset)) +
    geom_hline(yintercept = 0, linetype = 2, color = "gray70") +
    geom_pointrange(aes(ymin = lower, ymax = upper), size = 0.7, fatten = 0.5, position = position_dodge(width = 0.8)) +
    coord_flip() +
    scale_y_continuous(limits = c(-0.15, 0.3), breaks = c(-0.1, 0, 0.1, 0.2)) +
    facet_wrap(~plotorder, scales = "free_y", ncol = 1) +  # Uses plotorder to show 'position' label
    scale_colour_manual(values = c("gray80", "gray30"), labels = c("simple", "sophisticated")) +
    theme_fira() +
    labs(title = title_text, x = "", y = "") +
    theme(legend.title = element_blank(), legend.position = "bottom")
}

# Generate Plots Figure 8 - Left and Right Panels 
plot_mp_position(filter(mp_position_data, dataset %in% c("Low-Info Simple", "Low-Info Sophisticated")), 
                 "Low-Information environment")

plot_mp_position(filter(mp_position_data, dataset %in% c("High-Info Simple", "High-Info Sophisticated")), 
                 "High-Information environment")